In this notebook we use QGIS as a tool for interacting with data using NCI's web data services, and on the VDI filesystem.
We start by demonstrating one way of getting data into QGIS, and end with a shareable interactive 3D visualisation.
The program for this session is:
QGIS (Quantum GIS) is an application for geospatial data analysis and visualisation, focussed on 2D (raster and vector) datasets, and wrapped in a straightforward graphical user interface. It can also be used as a Python library, but for this session we use the GUI.
QGIS can be used on the VDI as a rapid-prototyping tool. If you've got all the infrastructure you need available in your project space and don't need it, that's great!
We hope to demonstrate some ways which QGIS can help you quickly develop shareable, multi-source analyses quickly and simply. We also revise some web service principles.
To answer our question we need some data on topography, vegetation and where block boundaries are.
From NCI's collections we can use:
NCI doesn't hold cadastral data, so we'll get those from the local land administrator:
Note - here, sections are collections of 50 or so individual house blocks. Our elevation and tree cover data are too coarse for individual house blocks - and the vector data file for blocks is massive!
For topography we need a terrain model in a raster format, e.g. GeoTIFF, which covers Canberra.
Shuttle Radar Topography Mission (SRTM) data are a good source of reasonably detailed elevation data held at NCI as NetCDF tiles.
We can head to the NCI Elevation collection here:
http://dapds00.nci.org.au/thredds/catalog/rr1/Elevation/catalog.html
...and look for SRTM 1 second elevation as NetCDF. Browse to the NetCDF folder, click through to:
http://dapds00.nci.org.au/thredds/catalog/rr1/Elevation/NetCDF/1secSRTM_DEMs_v1.0/DEM/catalog.html
A bunch of 1 degree tiles are shown here - but we will collect a region we want from the national mosaic:
Click on the 'WCS' link to see the WCS getCapabilities statement - which describes the data you can obtain here. We need to know the name for the coverage we need - look for the 'name' tag. With that,and using our WCS knowledge, we can request just the part of the data we need and acquire a GeoTIFF image. Let's break a WCS request down into parts we need:
To build a WCS request we concatenate the data path and dataset name, put a question mark after the dataset name, then add the rest of the labels describing the thing we want afterward, in any order, separated by ampersands:
Enter this link into a web browser to obtain a GeoTIFF DEM in your default download location (check in 'downloads'). Rename it to something memorable if you wish.
Note the use of GeoTIFF_Float - using only GeoTIFF is possible, but gives you an image with pixel values scaled to a colour range,m not a data range
Now go to QGIS and import the GeoTIFF as a raster layer:
We repeat this process with data on photosynthetically active vegetation, which comes from the ANU Water and Landscape Dynamics collection:
http://dapds00.nci.org.au/thredds/catalog/ub8/au/FractCov/PV/catalog.html
We'll look at 2015:
Using WCS again - click on the WCS link, and look for a
...which are assembled into a WCS request:
Again, using this link in a web browser will result in a GeoTIFF file arriving at your default download location. Load it into QGIS:
So far we have a bunch of grey things - let's make them colourful!
We don't need to worry about the elevation data, but we do need a useful colour scheme for vegetation data and our cadastral data. Double click your vegetation layer (in the layers panel) to bring up it's properties dialogue. From there:
(image above missing 'classify' step)
Steaming ahead? Add a graticule to your QGIS map layer...
How could you automagically save WCS data with sane, memorable names?
Why not just use the QGIS OWS browser to grab these data?
If you're keen, do some exploring - we don't know *everything* about QGIS - surprise us!
Why not get these data from the filesystem?
Great question! The main reason is that we're demonstrating QGIS as a way to fuse data from manifold sources. Over today and tomorrow you'll see a log examples of how to collect data from the filesystem in manageable chunks - which you could then analyse/visualise using the methods shown here. Demonstrating web coverage services on the VDI shows how you can pull data from many external sources to help interpret your work.
ACT publish a lot of spatial data on their ACTMAPi interface. To shortcut, here are the cadastral section data:
http://actmapi.actgov.opendata.arcgis.com/datasets/eedcda7873934e789093b093521b0299_3
You can download a shapefile and import it to QGIS, but let's show a feature you might like to use: adding vector data as a service. At ACTMAPi, click the 'API' menu box, and copy the URL out of the GeoJSON tetxtbox:
http://actmapi.actgov.opendata.arcgis.com/datasets/eedcda7873934e789093b093521b0299_3.geojson
In QGIS, Click 'add new vector layer', select the 'protocol' radio button, and past the URL in. The layer will load - but like raster data, it can't be used for analysis. Save it as a Geopackage and reload the layer (delete the GeoJSON layer).
Why Geopackage? Why not a shapefile?
Hint: http://www.geopackage.org
To procees we need to install two QGIS plugins - Zonal statistics and Qgis2threejs:
Do the same for a plugin named 'Qgis2threejs'.
We will try to show:
We start with Zonal statistics. This collects data from raster layers using polygons in a vector layer and computes basic statistics for each polygon.
As a proxy for section hilliness, we'll use the standard deviation of elevation in each block polygon.
Head to raster -> zonal stats to open the plugin.
In the zonal statistics plugin dialogue, choose the DEM as the raster layer, and use band 1. Then choose your ACT blocks vector layer. In the statistics to calculate, pick an appropriate set - but include standard deviation - this is our roughness proxy. Add a meaningful prefix to the statistics (e.g. 'dem_'), so you can find them when you need to use them.
QGIS will spend some time calculating stats for each block - grab a cup of tea!
For vegetation cover, each grid cell/pixel shows a percentage of photosynthetically active vegetation. As a quick measure we could sum all the pixel values inside a block to get an idea of it's vegetation coverage relative to other blocks.
In the zonal statistics plugin dialogue, choose the DEM as the raster layer, and use band 1. Then choose your ACT sections vector layer. In the statistics to calculate, the values we need are preselected - we'll use the mean value for each section. Again, use a meaningful prefix (for example, 'veg_')
We avoided styling our vector layer earlier, but now it's time - since we want to visualise the vegetation coverage in each section.
Double click your act_sections layer to open it's properties dialogue, and head to the style tab. Here, we want to apply a good looking colour scale based on mean vegetation cover (the data we just generated). Set up as follows:
...which results in the following map:
So far we can visualise the vegetation cover of sections in the ACT - but how can we relate that quickly and easily to section hillines? And how could we make a way to show the result in a convincing way?
is a graphical GIS the only way to collect zonal statistics using a vector layer to segment raster data?
So far we've imported three different datasets into QGIS and created some new attributes on our vector dataset. How
In this scheme, if hillier blocks are more vegetated, dark green blocks will visualise as taller columns. If hiller blocks are less vegetated, light green blocks will visualise as taller columns. Lets test it out!
Here we use the second plugin - Qgis2threejs. This renders the current screen to a WebGL map in a browser using three.js - with some neat data visualisation features.
So far we've worked in a WGS84 (EPSG:4326) coordinate system - but in order to render our map in three.js, we need to project our data into something which has units of metres, not degrees. Let's choose GDA94/MGA55 (EPSG:28355):
You'll see that everything has warped a touch, and your CRS panel (lower right) reflects your choice. Proj.4 handles all the rest for you!
Qgis2threejs attempts to render the whole map window as webGL. We want to limit or map to the extents of our DEM - so we can either zoom the map window in so that our region of interes occupies the whole window, or set a clipping polygon in a new vector layer. Here, we zoom in so that our region of interest fills the map window:
Head to web -> Qgis2threejs to open the plugin dialogue.
(some screebgrabs showing qgis2threejs setup)
In [6]:
###ignore this block of code - it is required only to show the map in iPython - you won't need it!
from IPython.core.display import display, HTML
display(HTML('<iframe width="800" height="600" frameborder="1" scrolling="no" src="./qgis2threejs/veg_by_hilliness.html"></iframe>'))
The folder you just created using QGIS2threeJS can be dropped into any web server - to share with collaborators.
You can't serve data directly from the VDI - but you can copy your results to some web host (eg github.io). Here's an example:
https://adamsteer.github.io/nci_samples/qgis2threejs/veg_by_hilliness.html
What are some other ways to share QGIS projects?
One option is to use curl on the command line:
curl -o SRTM_dem_139_36.tiff 'http://dapds00.nci.org.au/thredds/wcs/ub8/au/FractCov/PV/FractCover.V3_0_1.2015.aust.005.PV.nc?service=WCS&version=1.0.0&Request=GetCoverage&Coverage=PV&bbox=149.0,-36,149.9,-35&format=GeoTIFF_Float'
curl -o 2015_treecover_139_36.tiff 'http://dapds00.nci.org.au/thredds/wcs/ub8/au/FractCov/PV/FractCover.V3_0_1.2015.aust.005.PV.nc?service=WCS&version=1.0.0&Request=GetCoverage&Coverage=PV&bbox=149.0,-36,149.9,-35&format=GeoTIFF_Float'
...are there others? What was yours?
You could also try rasterstats: https://github.com/perrygeo/python-rasterstats - any other suggestions?
For simple projects with CCBY4 data, the QGIS cloud is one way of sharing your results. Another is a Jupyter notebook hosted on github as a gist, a notebook, or as a github.io page. What was your approach?
Loading OWS data is also possible in hosted 3D visualisers, for example NationalMap. Here's a sample of some XXXX data from NCI visualised using the NationalMap interface: (url) . What else do you use?
In [ ]: